function [t_stable_outside_SM,S_stable_outside_SM,...
    t_stable_inside_SM,S_stable_inside_SM,...
    t_unstable_inside_SM,S_unstable_inside_SM,...
    t_unstable_outside_SM,S_unstable_outside_SM] = PropagateManifold_L2(Xo,tau,u,numpoints,integrationTime)
epsilon=2*10^(-6);

T = tau*2;
%[] Single orbital period

to = linspace(0,T,numpoints);
%[]

options=odeset('RelTol',1e-14,'AbsTol',1e-22);
%[] Increase the Ode 45 tol

[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
t = t';
%[] Integrate states.

Monodromy = State2STM(t,S,u);
%[] Monodromy matrix

[eig_vec,eig_val] = eigs(Monodromy);
%[] Find eigen values and vectors

eig_val = diag(eig_val);
%[] Get just the value.

[~,ind_stable] = min(eig_val);
[~,ind_unstable] = max(eig_val);
%[] Determine stable and unstable index

nu_stable = eig_vec(:,ind_stable);
nu_unstable = eig_vec(:,ind_unstable);
%[] Stable and unstable direction perturbation

InitialConditions_stable_outside = zeros(6,length(t));
InitialConditions_stable_inside = zeros(6,length(t));
InitialConditions_unstable_outside = zeros(6,length(t));
InitialConditions_unstable_inside = zeros(6,length(t));

%[] Pre allocate initial conditions matrix

% Loop for creating initial conditions for manifold trajectories.
for ii = 1:(length(t)-1)
    
    t_STM = [t(1),t(ii+1)];
    %[] Determine the state transition matrix at desired time T
    
    Pert_unstable = State2STM(t_STM,Xo,u)*nu_unstable;
    Pert_stable = State2STM(t_STM,Xo,u)*nu_stable;
    %[] calculate the perturbance needed for unstable and stable manifold
    %trajectories
    
    Norm_Pert_unstable = Pert_unstable / norm(Pert_unstable);
    Norm_Pert_stable = Pert_stable / norm(Pert_stable);
    %[] Normalize the perturbing force.
    
    if Norm_Pert_stable(1)<0 % If x perturbation direction is negative
        
        InitialConditions_stable_outside(:,ii)      = S(:,ii) - epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) + epsilon * Norm_Pert_stable;
        disp('a')
    else
        InitialConditions_stable_outside(:,ii)      = S(:,ii) + epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) - epsilon * Norm_Pert_stable;
        disp('b')
    end
    
    if Norm_Pert_unstable(1)<0 % If x perturbation direction is negative
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) - epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) + epsilon * Norm_Pert_unstable;
        disp('a')
    else
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) + epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) - epsilon * Norm_Pert_unstable;
        disp('b')
    end
    %[] Determine the initial states for stable manifold's initial states.
    
end
T_manifold_stable_outside = [integrationTime,0]; % Integrating backward in time
%[] Integration time for stable manifold coming from outside

T_manifold_unstable_inside = [0,integrationTime];% Integrating forward in time
options_SecondBodyCrossing=odeset('RelTol',1e-14,'AbsTol',1e-22,'Events',@(t,S)HaltProp(t,S,u));

% Determine limit range

%Loop for propagating the calculated initial conditions
parfor ii = 1:length(InitialConditions_stable_outside)-1
    
    [t_stable_outside_SM{ii},S_stable_outside_SM{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
        InitialConditions_stable_outside(:,ii),options);
    S_stable_outside_SM{ii} = S_stable_outside_SM{ii}';
    for j = 1:length(S_stable_outside_SM{ii})
        x = S_stable_outside_SM{ii}(1,j);
        y = S_stable_outside_SM{ii}(2,j);
        theta = atan2d(y,x);
        projectedValue = Rotation([0;0;1],-theta,'Degrees')*S_stable_outside_SM{ii}(1:3,j);
        if projectedValue < (1-u)
            S_stable_outside_SM{ii} = NaN;
            break
        end
    end
    
    %[] Integrate the stable manifold trajectory coming in from outside.
    %
    [t_stable_inside_SM{ii},S_stable_inside_SM{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
        InitialConditions_stable_inside(:,ii),options);
    S_stable_inside_SM{ii} = S_stable_inside_SM{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
    %
    [t_unstable_inside_SM{ii},S_unstable_inside_SM{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
        InitialConditions_unstable_inside(:,ii),options);
    S_unstable_inside_SM{ii} = S_unstable_inside_SM{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
    
    [t_unstable_outside_SM{ii},S_unstable_outside_SM{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
        InitialConditions_unstable_outside(:,ii),options);
    S_unstable_outside_SM{ii} = S_unstable_outside_SM{ii}';
    for j = 1:length(S_unstable_outside_SM{ii})
        x = S_unstable_outside_SM{ii}(1,j);
        y = S_unstable_outside_SM{ii}(2,j);
        theta = atan2d(y,x);
        projectedValue = Rotation([0;0;1],-theta,'Degrees')*S_unstable_outside_SM{ii}(1:3,j);
        if projectedValue < (1-u)
            S_unstable_outside_SM{ii} = NaN;
            break
        end
    end
    %[] Integrate the stable manifold trajectory coming in from inside.
    
end
end